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Introduction 

The  determination  of  the  effect  of  surface  reflection  and  other  secondary  arrivals  on  the 
estimated  yields  of  nuclear  explosions  has  been  a  major  problem  in  nuclear  monitoring  work. 
The  major  philosophical  issues  in  such  work  are  the  following: 

a)  Are  the  source  functions  simple  superpositions  of  P  and  pP  pulses  or  are  more  com¬ 
plex  source  time  histories  appropriate? 

b)  Can  we  predict  the  effect  of  the  secondary  arrivals  on  the  yields  estimated  from  time 
domain  measurements  on  the  first  cycle  of  the  P  wave  signal? 

c)  What  is  the  effect  of  these  secondary  arrivals  on  frequency  domain  measures  of 
yield? 

d)  What  are  the  limitations  on  the  resolution  of  secondary  arrivals  imposed  by  Q  and 
the  noise? 

e)  And,  in  general,  what  processes,  e.g.  Q,  scattering,  near  source  and  receiver  wave 
conversions,  spall,  Rayleigh-to-P  conversions,  etc.,  determine  the  waveforms  and 
spectra  of  teleseismic  P  waves? 

These  problems  are  not  trivial.  In  the  case  of  a  simple  pP  reflection  we  understand  how 
it  biases  the  various  time  domain  measurements  and  we  are  thus  able  to  make  corrections  for  it 
if  the  amplitude,  arrival  time  and  spectral  characteristics  of  the  pP  are  known.  If,  on  the  other 
extreme,  the  P  wave  consists  of  various  multipathed  and  focused  arrivals  in  conjunction  with 
complex  source  functions  consisting  of  many  pulses  of  unknown  origin,  then  there  is  no  physi¬ 
cal  model  on  which  we  could  base  any  methods  to  "correct"  our  readings,  since  it  becomes 
doubtful,  not  understanding  the  nature  of  arrivals,  that  even  the  first  arrival  amplitude  is  uncon¬ 
taminated  and  well  understood.  The  pP  yield  correction  is  based  on  the  belief  that  there  are 
only  two  pulses  present,  the  amplitude  of  the  first  of  which  is  a  measure  of  direct  radiation 
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from  an  explosive  source.  If  there  are  other  pulses  present  then  we  still  have  to  assume  that 
the  first,  upward  motion  is  due  to  an  uncontaminated  direct  arrival  from  the  source. 

Another  important  issue  is  the  basic  resolution  attainable  in  determining  the  time  delay  of 
pP,  provided  it  is  present.  For  band-limited  signals  and  small  pP  delays  we  may  not  be  able  to 
resolve  this  time  accurately  enough  to  get  a  meaningful  correction,  even  if  the  simple  P+pP 
model  is  valid.  For  such  signals  it  may  be  possible  to  obtain  spectacular  waveform  fits,  but 
these  do  not  prove  that  the  model  is  valid.  Most  of  these  questions,  the  validity  of  the  models, 
the  resolution  of  the  data,  and  the  nature  of  multiple  pulses,  have  to  be  decided  on  the  basis  of 
analyses  of  data  rather  than  by  intuitive  reasoning  or  fitting  the  data  to  unverified  theoretical 
models. 

In  this  report,  we  discuss  work  done  to  help  resolve  the  questions  discussed  above  about 
the  nature  of  the  first  several  seconds  of  the  direct  P  waveform  from  explosion  sources.  First, 
we  discuss  previous  work  on  this  problem  and  pitfalls  associated  with  attempting  to  identify 
the  pP  phase  on  short  period  records  of  explosions.  Subsequently  we  shall  briefly  present 
results  from  spectral  ratio  techniques  for  estimating  the  pP-P  delay  times  and  amplitudes  and 
from  stacked  deconvolutions  which  enhance  the  pP  phase  on  seismic  records.  The  rest  of  this 
report  is  a  description  of  the  results  of  using  a  maximum  likelihood  multichannel  deconvolution 
technique  for  estimating  time  domain  spike  series  for  both  the  source  and  site  terms.  Four 
appendixes  at  the  back  of  the  report  contain  further  examples  of  the  data  and  results  of  the 
work.  Finally,  we  present  suggestions  for  future  work. 

We  find  that  explosion  source  time  functions,  while  they  generally  appear  to  contain  a 
phase  that  is  readily  ascribable  to  pP,  also  contain  a  number  of  other  arrivals.  Both  source  and 
site  terms  contribute  to  the  presence  of  P  coda.  The  deconvolution  of  a  presumed  cratering 
event  shows  a  more  complex  initial  P  waveform  than  that  obtained  from  buried  explosions, 
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without  a  clearly  discernible  pP.  From  these  deconvolutions,  we  estimate  the  cratering- 
noncratering  ir^  bias  to  be  0.09  to  0.22  magnitude  units. 

Background 

A  variety  of  methods  have  been  proposed  for  estimating  pP  times  and  amplitudes.  Spec¬ 
tral  methods  (Cohen  1970,  1975;  Shumway  and  Blandford  1978)  attempt  to  make  use  of  the 
spectral  modulation  due  to  pP  and  are  successful  only  if  no  other  secondary  arrivals  are 
present.  In  the  presence  of  numerous  secondary  arrivals  the  spectra  are  modulated  according  to 
combinations  of  all  of  the  time  differentials  of  the  various  arrivals,  thus  resulting  in  a  complex 
modulation  pattern  that  cannot  be  inverted  easily.  In  the  case  of  complex  site  related  modula¬ 
tion  of  P  wave  spectra,  the  spectral  minima  due  to  any,  even  solitary,  pP  arrival  besides  the 
direct  P,  may  be  completely  disguised,  thus  making  the  interpretation  of  spectral  minima 
difficult.  Moreover,  even  in  the  absence  of  site  related  modulations  it  is  impossible  to  recon¬ 
struct  the  source  related  spike  sequences  from  power  spectra,  since  these  methods  discard  the 
phase  information.  The  same  power  spectrum  can  correspond  to  an  infinite  set  of  spike 
sequences  and  these  spike  sequences  can  be  derived  from  each  other  by  applying  arbitrary 
"dispersive"  filters  to  them  (Robinson  and  Treitel  1980).  In  Figures  1  and  2  we  illustrate  sets 
of  time  functions  that  have  exactly  the  same  power  spectra.  Several  of  these  could  be  inter¬ 
preted  as  complex  npF'  and  "spall"  arrivals  following  a  P.  As  soon  as  we  allow  even  a  dis¬ 
torted  pP  we  have  nonuniqueness  since  a  combination  of  a  "P”  and  a  "low  frequency  pP"  can 
have  spectra  not  different  from  that  of  a  single  delta  function  (Figure  1). 

Time  domain  methods  include  various  deconvolution  schemes  (Marshall  1972;  Bakun  and 
Johnson  1973;  Kemerait  and  Sutton  1982;  Israelson  1983)  which  generally  preserve  the  phase 
information.  Least  squares  (12  norm)  deconvolution  models  work  well,  but  those  that  use  the  lj 
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norm  (Levy  and  Fullagar  1981;  Mellman  et  a]  1984)  were  shown  by  Jurkevics  and  Wiggins 
(1984)  to  result  in  erroneous  deconvolutions  as  judged  by  synthetic  simulations.  Another 
method  that  has  been  tried  for  the  determination  of  pP  times  is  the  direct  modeling  of  P 
waveforms.  It  has  been  shown  (Cormier  1982)  that  such  waveform  fits  are  highly  nonunique 
because  t*,  pP  times  and  amplitudes,  and  the  various  explosion  source  parameters  trade  off 
against  each  other.  A  different  approach  that  has  been  proposed  involves  cross-equalization  (or 
as  the  authors  call  it,  "intercorrelation”)  of  waveforms  (Mellman  and  Kaufman  1981;  Lay  et  al 
1984a).  This  method  produces  "pseudo-pP"  delay  times  and  amplitudes,  i.e.  the  best  fits  to 
simple  P+pP  models.  Since  the  results  are  not  deconvolutions  in  the  true  sense  they  do  not 
prove  the  existence  or  correctness  of  such  models.  The  problem  with  such  time  domain 
methods  is  that  there  is  no  established  quantitative  relationship  between  the  quality  of  fit, 
where  the  fit  is  commonly  based  on  the  use  of  some  kind  of  a  norm  (Lay  et  al  1984a;  Mellman 
et  al  1984),  and  the  level  of  confidence  attached  to  a  hypothesis  that  the  P+pP  model  is  valid. 
This  can  be  illustrated  by  synthetically  simulating  the  case  in  which  many  secondary  arrivals 
are  present.  We  have  constructed  "receiver  functions"  by  appending  an  impulse  with  a  gradu¬ 
ally  decaying  Gaussian  noise  sequence.  For  the  source  pulse  we  have  used  a  narrow  band 
pulse  of  the  type  commonly  generated  by  seismic  sources  at  teleseismic  distances  as  seen 
through  WWSSN  type  instruments  and  then  generated  a  "source  spike  series"  by  following  a 
unit  impulse  with  a  burst  of  Gaussian  random  noise  of  unit  maximum  amplitude  with  a  dura¬ 
tion  of  1.25  seconds.  These  are  certainly  more  complex  than  a  P+pP  sequence.  Despite  this, 
applying  the  "intercorrelation"  method  we  were  able  to  achieve  correlation  coefficients  of  the 
same  order  as  those  seen  in  the  study  of  Lay  et  al  (1984a).  Figures  3  and  4  show  four  exam¬ 
ples  of  the  input  traces,  the  traces  after  intercorrelation,  and  the  resultant  correlation  coefficient. 
Thus  the  "intercorrelations"  with  real  data  prove  that  the  P+pP  models  are  too  simplistic,  since 
otherwise,  instead  of  values  near  0.6  to  0.9,  we  should  be  able  to  obtain,  routinely,  correlation 
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Figure  3.  Synthetic  inputs  (on  left)  and  their  "intercorrelations"  (on  right)  with  the  resulting 
correlation  coefficients. 
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coefficients  near  0.95.  It  appears  that  with  narrow  band  data  it  is  possible  to  convert  practi¬ 
cally  anything  into  anything  as  long  as  the  spike  sequences  of  the  various  sources  have  short 
time  durations  (which  seems  to  be  the  case).  Thus,  the  intercorrelation  method  may  sometimes 
yield  valid  results,  but  only  if  the  obtained  correlation  coefficients  are  very  high,  say  above 
0.95.  Most  of  the  time,  however  the  low  values  of  the  correlation  coefficients  obtained  after 
"intercorrelation"  indicate  that  the  underlying  model  is  wrong. 

As  the  last  item  on  our  list  of  deconvolution  methods,  we  mention  the  approach  sug¬ 
gested  by  Mellman  et  al  (1984)  that  consists  of  a  set  of  trial  deconvolutions  varying  the  P-pP 
time  lag.  It  is  claimed  that  the  changes  in  the  waveform  characteristics  and  amplitudes  in  the 
deconvolution  results  indicate  the  "best"  pP  delay.  Unfortunately,  this  is  not  the  case.  Since 
the  inverse  of  the  P+pP  pair  is  an  exponentially  decreasing  comb  filter  it  will  always  produce 
the  simplest  waveform  with  the  minimum  amplitude  near  half  the  dominant  period  of  the  origi¬ 
nal  wave  and  will  give  a  ringing,  large  amplitude  waveform  for  lag  times  close  to  the  period  of 
the  wave.  The  nature  of  the  output  is  thus  determined  by  the  dominant  period  and  has  nothing 
to  do  with  pP  characteristics. 

Basically,  the  same  kinds  of  problems  affect  all  methods  proposed  for  the  estimation  of 
the  pP  times  of  nuclear  explosions.  Generally  speaking,  it  may  be  stated  that  all  proposed 
techniques  seem  to  work  quite  well  with  synthetic  data  that  incorporate  the  assumptions  about 
the  makeup  of  the  P  wavetrain  that  were  used  in  the  design  of  the  deconvolution  method  itself. 
Most  of  the  new  methods  for  determining  pP  parameters  have  been  carefully  checked  to  verify 
that  the  method  works  for  the  signal  model  it  was  based  on.  The  same  care  was  not  applied, 
generally,  to  verify  from  the  data  that  the  underlying  models  are  valid.  This  explains  why  the 
results  obtained  from  real  data  are  so  ambiguous.  The  advantage  of  the  maximum  likelihood 
estimation  method  discussed  later  in  this  study  is  that  it  assumes  little  about  the  nature  of  the 
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source  and  site  factors,  and  it  is  based  mostly  on  the  factorable  nature  of  P  wave  spectra  which 
is  verified  by  the  successful  reconstruction  of  the  original  data.  We  must  caution  the  reader 
also  that  the  successful  reconstructions  (resulting  in  high  values  of  measures  of  coherence 
between  synthetic  and  real  data)  or  cross-equalizations  (Lay  et  al  1984a)  of  narrow-band  body 
wave  data  do  not  by  themselves  prove  the  correctness  or  uniqueness  of  the  pP,  spall,  and  other 
parameters  derived  by  such  methods.  Generally,  waveform  matches  are  nonunique  and  are  not 
valid  substitutes  for  true  deconvolutions  in  which  the  multiple  arrivals,  if  any,  are  separated 
and  become  identifiable  by  the  enhancement  of  the  high  frequency  content.  Only  deconvolved 
records  overcome  the  weakness  in  human  perception  that  makes  us  unable  to  distinguish  multi¬ 
ple  arrivals  convolved  with  filters  that  severely  limit  the  bandwidth.  In  addition,  statistical 
techniques  are  needed  to  gauge  the  resolution  available  with  the  given  data. 

We  conclude,  therefore,  that  our  preconceptions  about  the  structure  of  teleseismic  body 
waves  are  more  likely  to  be  at  fault  than  any  techniques  that  were  proposed.  In  other  words,  if 
the  teleseismic  P  waves  were  simply  convolutions  of  a  P-pP  pair  of  pulses  with  some  random 
time  function  appropriate  to  o  seismic  path  the  problem  of  finding  unique  pP  estimates  would 
have  been  solved  long  ago.  The  continuing  "lack  of  success"  in  this  field  must  thus  be  due  to 
the  fact  that  even  "simple"  events  such  as  nuclear  explosions  radiate  quite  complex  wavetrains 
and  the  underlying  physical  processes  are  not  well  understood. 

Spectral  Ratios 

In  this  section,  we  briefly  discuss  the  results  of  using  a  simple  spectral  ratio  method  to 
put  constraints  on  the  pP  amplitude  and  a/nval  time.  Spectral  ratios  between  two  explosions 
from  the  same  source  region  were  taken  at  various  stations  of  NORSAR.  Thus,  since  the 
events  were  from  the  same  source  region,  near-source  site  effects  should  nearly  cancel  out,  as 
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should  receiver  effects  since  both  events  in  each  pair  arrive  on  nearly  the  same  azimuth.  The 
sequence  of  peaks  and  nulls  that  make  up  the  spectra  should  then  simply  be  due  to  the  interfer¬ 
ence  of  the  different  pP  amplitudes  and  delay  times  of  the  two  events. 

Representative  spectral  ratios  between  the  same  pair  of  events  at  different  NORSAR 
receivers  are  shown  in  Figures  A1  to  A4  of  Appendix  A.  Synthetic  spectral  ratios  are  also 
shown  on  each  figure.  These  were  made  assuming  a  von  Seggem  and  Blandford  (1972)  source 
time  function  for  each  event,  and  using  yield,  pP  amplitude,  and  pP  time  estimates  from  Lay 
et  al  (1984a)  as  shown  in  Table  A1  in  Appendix  A.  For  each  pair  of  events,  the  spectral  ratios 
at  the  different  receivers  have  many  similarities  in  even  small  details  of  these  spectral  ratios  at 
the  various  receivers.  It  is  clear  that  the  synthetic  spectral  ratios  derived  from  the  simple  P+pP 
model  are  generally  not  consistent  with  the  observed  ratios.  This  cannot  be  attributed  merely 
to  incorrect  estimates  of  pP  delay  times  and  amplitudes,  but  a  more  likely  reason  is  the  pres¬ 
ence  of  additional  arrivals  in  the  initial  part  of  the  P  wave  record  besides  possible  P  and  pP. 
The  existence  of  complex  multiple  arrivals  originating  from  the  source  is  also  confirmed  by 
our  deconvolution  results  presented  below. 

The  presence  of  additional  complications  in  the  source  time  function  is  also  suggested  by 
the  spectra  themselves.  Few  spectra  from  nuclear  explosions,  with  some  notable  exceptions 
such  as  the  very  deep  Amchitka  shots  and  the  Bukara  shot  (Marshall  1972),  show  clear  spec¬ 
tral  nulls  at  the  times  predicted  from  the  interference  of  P  and  pP.  Also,  spectra  of  synthetic 
wavelets  generated  for  propagation  through  random  media  by  the  finite  difference  method  show 
only  hints  of  spectral  nulls,  while  spectra  of  synthetics  generated  by  propagation  through  sim¬ 
ple,  plane  layered  media  show  distinct  spectral  nulls  (McLaughlin  et  al  1986).  An  example  of 
such  results  is  shown  in  Figure  5,  which  compares  the  expected  teleseismic  P  waveforms  from 
explosion  sources  in  a  two-dimensional  model  of  Yucca  Flats  with  those  computed  for  a  simi- 
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Figure  5.  (a)  West-to-east  model  for  geologic  structure  across  Yucca  Valley  used  in  finite 
difference  simulations.  The  model  is  shown  with  5-to-l  vertical  exaggeration,  and  each  region 
is  labeled  with  its  P-wave  velocity.  P-wave  velocities  of  1.34,  2.14,  3.00,  and  4.57  km/sec  are 
indicated  for  the  geologic  units  of  alluvium,  unsaturated  tuff,  saturated  tuff,  and  Paleozoic  car¬ 
bonates,  respectively.  The  source  locations  are  indicated  by  numbered  dots  across  the  middle 
of  the  figure,  (b)  Synthetics  for  source  locations  4,  5,  and  6  at  take-off  angle  of  15  degrees 
(solid  lines)  compared  to  synthetics  for  a  plane-layered  model  (dashed  lines),  (c)  Spectra  of 
the  synthetics  for  source  location  5  of  the  Yucca  Valley  model  (left)  and  for  the  plane-layered 
model  (right)  showing  the  much  more  distinct  spectral  nulls  in  the  spectra  of  the  synthetics 
from  the  flat-layered  model  (after  McLaughlin  et  al  1986). 
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lar  but  flat  layered  structure.  After  the  initial  1.5  seconds,  the  waveforms  from  the  two  models 
are  very  different.  This  example  illustrates  that  complex  multiple  arrivals,  not  simple  combina¬ 
tions  of  P  and  pP,  should  be  generally  observed  from  many  laterally  varying  geological  struc¬ 
tures  containing  sources.  The  physical  nature  of  the  complex  later  arrivals  in  this  case  is  prob¬ 
ably  conversion  of  Rayleigh  waves  to  P  at  the  sides  of  the  valley. 

Complications  are  also  observed  in  the  receiver  terms.  Figures  A5  to  A8  in  Appendix  A 
show  intersensor  spectral  ratios  at  NORSAR  for  Nevada  Test  Site  (NTS)  nuclear  explosions; 
use  of  ratios  cancels  out  the  source  spectra.  While  the  shapes  of  spectral  ratios  below  4  Hz  are 
similar  for  the  same  pairs  of  sensors,  the  strong  variations  in  these  ratios  with  frequency  indi¬ 
cates  that  details  in  the  spectral  shapes  at  any  given  site  are  unstable  for  arrivals  from  a  limited 
source  region.  This  result  is  not  surprising,  since  waves  propagated  through  media  with  ran¬ 
dom  velocity  variations  commonly  reveal  variations  in  details  of  their  spectra  while  the  gross 
shapes  are  relatively  constant  after  propagating  roughly  the  same  distance  (Frankel  and  Clayton 
1984).  Moreover,  even  if  one  assumes  a  flat  layered  medium  for  a  crustal  structure,  it  would 
change  the  details  in  the  spectral  shapes,  while  leaving  the  overall  slope  relatively  constant. 
Thus,  there  are  too  many  variations  in  the  source  and  receiver  information  derived  from  the 
spectral  ratio  methods  used  here  to  provide  significant  estimates  of  the  amplitude  and  delay 
time  of  pP  or  to  provide  information  on  other  phases  which  are  part  of  the  initial  P  wavetrain. 

Stacked  Deconvolutions 

Another  simple  method  for  the  estimation  of  pP  parameters  is  by  deconvolving  the 
records  assuming  that  the  source  pulse  is  known.  We  assume  that  the  source  pulse  is  made  up 
of  the  explosion  far-fleld  pulse,  a  constant  Q  operator,  and  the  instrument  impulse  response. 
Although  several  models  were  proposed  for  the  explosion  far  field  pulse  and  its  scaling  with 
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respect  to  yield  and  depth  (Mueller  and  Murphy  1971;  von  Seggern  and  Blandford  1972,  Lay 
et  al  1984b),  the  differences  between  them  do  not  influence  the  practical  results  much.  To 
illustrate  this,  in  Figure  6  we  show  some  deconvolutions  of  wavelets,  where  the  wavelets  were 
derived  from  convolving  a  short  period  instrument  response  and  an  attenuation  operator  with 
source  pulses  of  the  von  Seggem-Blandford  and  Mueller-Murphy  types.  The  deconvolutions 
used  two-sided  Wiener  filters  designed  from  the  autocorrelation  functions  of  the  wavelets  to 
which  a  50%  delta  function  was  added  to  simulate  white  noise.  The  effect  of  this  adjustment 
is  to  make  the  deconvolution  less  precise,  thus  avoiding  swamping  the  results  with  amplified 
high  frequency  noise.  This  has  the  same  effect  as  the  0  parameter  in  the  maximum  likelihood 
deconvolutions  shown  below.  Each  wavelet  was  deconvolved  separately  with  both  source 
pulses.  The  only  effect  of  using  the  wrong  type  of  source  pulse  in  the  deconvolutions  is  some 
asymmetry  in  the  deconvolved  wavelets  which  would  hardly  cause  problems  in  any  practical 
work  in  separating  pulses.  Therefore,  the  most  important  contributions  to  the  shape  of  the 
"wavelet"  are  the  instrument  response  and  the  Q  operator.  While  the  explosion  source  model 
includes  a  spectral  scaling  feature,  which  helps  us  to  deal  with  explosions  of  vastly  different 
sizes,  we  could  have  used  an  impulse  instead  and  the  results  would  not  have  been  much 
different.  Thus  it  appears  that  the  past  lack  of  success  in  deconvolutions  (or  rather  in  finding  a 
simple  pP)  of  real  data  could  hardly  be  due  to  uncertainties  in  the  explosion  source  models. 
Even  if  we  had  used  the  model  by  Lay  et  al  (1984b),  a  model  that  in  our  opinion  is  invalid 
because  it  gives  the  wrong  high  frequency  spectral  falloff  rate  (Murphy  and  Day  1984),  it 
would  not  have  affected  the  results  much. 

The  time  domain  representation  of  the  attenuation  operator  also  depends  somewhat  on 
whether  Q  is  frequency  dependent  or  not  in  the  short  period  band,  but  refinements  such  as 
those  to  our  models  must  wait  until  more  is  known  about  the  frequency  dependence  of  Q.  The 
apparent  t*  along  most  seismic  paths  should  be  estimable  from  the  gross  shapes  of  spectra  with 
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Figure  6.  Deconvolutions  of  theoretical  source  time  functions.  The  original  source  time  func¬ 
tions  on  the  left  were  made  using  the  von  Seggem  and  Blandford  (1972)  model  (vSB)  while 
those  on  the  right  were  made  using  the  Mueller  and  Murphy  (1971)  model;  all  four  source  time 
functions  were  deconvolved  with  the  vSB  model  using  a  Wiener  filter.  The  upper  two  source 
functions  were  not  prewhitened;  the  lower  two  were  prewhitened  by  increasing  the  zero  lag 
autocorrelation  value  by  50%.  The  only  difference  in  the  deconvolutions  due  to  the  different 
source  functions  is  a  slight  asymmetry. 
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an  accuracy  of  about  0.1  sec  (Der  et  al  1985). 

In  order  to  reduce  the  effects  of  site  related  distortion  of  the  pulses  we  follow  Marshall 
(1972)  in  stacking  the  deconvolutions  over  widely  distributed  sensors  of  large  arrays.  This 
way  we  obtain,  hopefully,  deconvolved  source  spike  sequences  in  which  the  site  effects  are 
averaged  out,  or  at  least,  reduced.  These  results  are  equivalent  to  the  outputs  of  the  first  itera¬ 
tion  in  the  multichannel  deconvolution  method  described  later  in  this  report.  As  we  shall  see 
later  the  results  are  quite  similar,  but  because  of  the  more  refined  nature  of  the  latter  we  place 
more  trust  in  the  results  of  the  iterative  method.  Nevertheless,  the  results  of  this  simpler 
approach  serve  as  good  preliminary  checks  of  the  multichannel  iterative  method. 

Results  of  the  stacked  deconvolutions  for  5  NTS  events  at  NORSAR  are  shown  in  Figure 
7.  Most  of  these  events  show  a  late  negative  pulse  near  1  second  after  the  first  arrival,  but  are 
more  complex  than  a  simple  P+pP  pair.  Stilton  shows  an  early  arrival  of  negative  polarity,  and 
Inlet  shows  a  very  early  arrival  of  negative  polarity. 

The  stacked  deconvolutions  of  four  Eastern  Kazakh  explosions  at  EKA  in  Figure  8  show 
an  early  (0.4  sec)  inverted  pulse  following  the  first  positive  pulse.  The  top  trace  appears  to  be 
different  from  the  rest  in  the  deconvolved  waveform;  this  is  an  assumed  cratering  event  that 
will  be  discussed  again  in  more  detail  below. 

The  stacked  deconvolutions  in  Figure  9  of  the  same  Kazakh  explosion  recorded  at  three 
AWRE  arrays  are  interesting  in  that  they  are  very  similar.  It  appears  that  this  event  was  fairly 
symmetric  azimuthally  and  was  probably  not  associated  with  significant  strain  release  or  other 
asymmetries  in  radiation,  at  least  in  the  short  period  band. 
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Figure  7.  Stacked  deconvolutions  of  five  NTS  events  recorded  at  NORSAR.  Estimates  of  the 
pP  delay  time  inferred  from  the  deconvolutions  are  shown  to  the  right  of  the  traces. 


15  Jan  1965 
Cratering  event 

T  =  .3  sec 
T  =  .7  sec 


27  April  1975 

T  =  .3  sec 
T  =  .6  sec 


Figure  8.  Stacked  deconvolutions  of  five  Eastern  Kazakh  explosions  recorded  at  EKA.  The 
top  event  is  a  presumed  cratering  explosion  and  has  a  noticeably  more  complex  deconvolved 
waveform.  Estimates  of  the  pP  delay  time  inferred  from  the  deconvolutions  are  shown  to  the 
right  of  the  traces. 


] 


Figure  9.  Stacked  deconvolutions  of  an  Eastern  Kazakh  explosion  (20  February  1975) 
recorded  at  the  AWRE  arrays  EKA,  GBA,  and  YKA.  The  similarity  of  the  results  for  the 
same  event  recorded  at  different  azimuths  suggests  that  the  source  was  fairly  symmetric.  Esti¬ 
mates  of  the  pP  delay  time  inferred  from  the  deconvolutions  are  shown  to  the  right  of  the 
traces. 
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Maximum  Likelihood  Multichannel  Deconvolutions 

Description  of  the  Deconvolution  Method 

A  major  advance  in  the  understanding  of  the  structure  of  teleseismic  P  waves  was  the 
finding  by  Filson  and  Frasier  (1972)  that  P  wave  spectra  from  a  set  of  events  observed  at  a 
narrow  range  of  azimuths  at  multiple  array  sites  may  be  factored  as 


Yjj((o)  =  Rj(e»  Sj(co)  ,  (1) 

where  Yjj(co)  are  the  P  wave  spectra,  Sj((0)  are  source  factors,  and  Ri(co)  are  factors  appropriate 
to  the  receiver  sites.  The  validity  of  this  equation  has  been  verified  for  a  number  of  seismic 
arrays  as  well  as  more  widely  separated  receivers  (Lundquist  et  al  1980)  and  the  idea  was  also 
exploited  for  deriving  "relative  receiver  functions",  which  are  time  domain  representations  of 
intersite  transfer  functions  (Lundquist  et  al  1980)  including  Q  effects. 

In  this  report  we  shall  talk  of  "receiver  factors"  and  "source  factors"  in  order  to  make  the 
distinction  from  the  related  concept  of  "relative  receiver  functions",  which  also  incorporate 
some  additional  assumptions  (Hart  et  al  1979),  such  as  the  maximum  "spikiness"  condition 
(Wiggins  1978),  which  are  not  required  by  the  physics  of  the  problem.  We  separate  the  Q 
operator,  the  instrument,  and  a  von  Seggem  and  Blandford  (1972)  source  wavelet  from  the 
"source  factor"  and  combine  these  in  an  assumed  known  source  wavelet,  Aj.  Thus, 
Sj((0)  =  Aj(co)Xj(o))  where  Aj(co)  is  the  spectrum  of  the  combined  source  wavelet,  instrument, 
and  t*.  This  way  the  "source  factor",  Xj(<a),  is  a  Fourier  transform  of  a  spike  sequence  associ¬ 
ated  with  the  radiation  from  a  source  and  the  "receiver  factor",  Rj(w),  is  the  Fourier  transform 
of  the  site  transfer  function.  Since  the  receiver  factors  strongly  depend  on  the  azimuths  and 
distances  of  the  sources,  they  may  be  assumed  to  be  constant  for  a  given  sensor  only  for  lim¬ 
ited  source  regions. 
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If  we  have  N  receivers  and  M  sources,  we  may  describe  (  N  x  M  )  P  wave  spectra  (or 
time  series)  in  terms  of  N  +  M  factors;  there  is  thus  a  considerable  degree  of  redundancy  in 
the  problems  with  N  and  M  sufficiently  large.  This  helps  us  to  separate  the  Rj(o>)  from  the 
Xj(to).  This  cannot  be  achieved  uniquely,  however,  unless  something  is  known  a  priori  about 
the  properties  of  each.  Otherwise,  common  spectral  factors  can  be  assigned  to  either  the 
source  or  receiver  factors,  or  canceling,  reciprocal  factors  may  be  assigned  to  both  the  receiver 
and  source  factors  as  indicated  below: 

Yij(co)  =[R;(co)  F(£0)1  [F-'tm)  X/0))]  (2) 


Since  all  P  wave  signals  are  embedded  in  the  noise  backgrounds,  it  is  desirable  that  the 
deconvolution-factoring  process  is  optimized  in  some  way  with  respect  to  the  noise  i.e.,  we  are 
trying  to  find  some  compromise  between  the  resolution  and  the  S/N  ratio.  Although  the  origi¬ 
nal  derivation  of  Shumway  (1984)  includes  weighting  with  respect  to  noise  spectra,  in  this 
work  we  are  not  striving  for  maximizing  the  S/N  ratio,  but  seek  some  compromise  between  the 
time  resolution  of  the  deconvolution  and  the  output  noise  level. 

The  detailed  derivation  of  the  multichannel  deconvolution  method  is  given  in  Der  et  al 
(1983)  and  Shumway  and  Der  (1985).  Only  the  key  formulas  are  repeated  here: 


X  Aj  R;  Yjj 

_ | _ 

X  I  Aj  |2  |  R;  |2  +  8j 


where  0j  =  ——  , 


(3) 


X  I  Aj  |2  I  Rj  I2  +  0j 


-i 


and 


(4) 
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where 


EAjXjY* 

_ j  _ 

D  M2  [|  Xj  i2  +  dj2] 

j 


(5) 


Yjj  data 

Aj  source  wavelet  spectrum 
Rj  receiver  spike  response 
Xj  source  spike  sequence 
0  adjustable  S/N  ratio 
Pv  noise  power  spectrum 
Pz  initial  signal  power 
Oj2  effective  noise  power  in  second  step 

and  an  overbar  denotes  the  complex  conjugate.  The  i’s  correspond  to  the  individual  receiver 
sites  and  the  j’s  to  the  individual  sources.  Equations  3  through  5  are  very  similar  to  those 
presented  by  Oldenburg  (1981)  and  Oldenburg  et  al  (1981),  and  are  essentially  multichannel 
versions  of  the  same. 

The  deconvolution  is  done  in  the  frequency  domain.  The  first  step  is  to  factor  out  the 
"known"  spectra  consisting  of  t*.  the  source  function,  and  the  instrument  response.  Since  we 
are  working  with  source  and  receiver  arrays,  we  assume  that  the  attenuation  parameter  t*P  is 
the  same  for  all  sensors  and  that  all  instruments  have  the  same  response.  At  the  start  of  the 
iteration  process  the  source  and  site  functions  are  each  initialized  to  an  impulse  in  the  time 
domain.  The  flow  diagram  of  the  iterations  is  depicted  in  Figure  10.  After  several  iterations, 
each  source  term  will  contain  the  terms  common  to  that  source  at  all  receivers  including  the 
source  pulse,  pP  and  spall  arrivals,  and  any  scattering  contributions  that  may  be  associated  with 
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Figure  10.  Flow  diagram  of  the  iterative  maximum  likelihood  multichannel  deconvolution 
algorithm.  In  the  center,  the  second  and  third  boxes  from  the  top  correspond  to  equations  (3) 
and  (5)  in  the  text,  respectively. 
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a  particular  source.  Each  site  term  will  contain  the  terms  common  to  that  site  for  all  sources, 
mostly  the  random  scattering  effects  attributable  to  the  structure  near  the  receiver.  The  results 
are  checked  by  reconstructing  each  input  trace  from  the  convolution  of  the  appropriate  source 
and  site  terms  with  the  wavelet  consisting  of  the  von  Seggem  and  Blandford  source  pulse,  t*, 
and  instrument  response.  The  differences  between  the  original  and  reconstructed  seismograms 
can  be  inferred  as  a  measure  of  the  scattering  unique  to  each  individual  path  i.e.,  energy  unac¬ 
counted  for  by  either  site  or  source  factors,  and  thus  missing  in  the  reconstructions. 

This  multichannel  deconvolution  approach  for  estimating  source  and  site  region  charac¬ 
teristics  has  several  advantages.  First,  the  method  utilizes  both  amplitude  and  phase  informa¬ 
tion  in  the  frequency  domain  calculations.  Furthermore,  no  a  priori  assumptions  need  to  be 
made  about  the  complexity  of  either  the  source  or  the  site  spike  series,  e.g.  the  presence  of  pP 
does  not  need  to  be  presupposed  by  providing  an  initial  guess  of  the  pP  delay  time  and  ampli¬ 
tude. 

These  iterative  formulas  yield  progressively  improved  estimates  of  the  source  and  receiver 
factors  and  their  equivalent  time  functions.  It  is  apparent  from  equation  (3)  that  the  procedure 
is  similar  to  the  iterative  beaming  method  (Blandford  et  al  1973)  in  which  the  previous  esti¬ 
mates  of  the  site  and  receiver  factors  are  deconvolved  in  each  step.  The  resolution  kernels 
characterizing  the  attainable  resolution  with  the  given  bandwidth  are  also  plotted  at  the  conclu¬ 
sion  of  the  iterative  process. 

Results  of  the  Data  Analyses 

Four  data  sets  have  been  analyzed  in  this  study:  five  nuclear  explosions  from  the  Nevada 
Test  Site  (NTS)  recorded  at  five  stations  of  NORSAR,  eight  explosions  from  the  East  Kazakh 
Test  Site  (KTS)  recorded  at  five  stations  of  NORSAR,  four  Astrakhan  explosions  recorded  at 
four  RSTN  stations,  and  5  KTS  explosions,  including  a  presumed  cratering  explosion,  recorded 
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at  twelve  stations  of  the  AWRE  array  EKA.  The  events  in  each  data  set  are  described  in 
Tables  1  to  4  and  the  parameters  used  in  the  deconvolutions  are  in  Table  5. 

The  following  approach  was  used  to  prepare  the  data.  All  traces  were  lined  up  in  time 
on  the  first  arrival.  Since  we  are  not  interested  in  the  absolute  amplitudes  in  the  context  of 
estimating  relative  pP  times  and  amplitudes,  and  because  v/e  want  to  give  roughly  equal  weight 
to  all  the  traces  to  ensure  effective  factoring,  we  multiplied  each  trace  with  two  factors,  one 
appropriate  to  the  event,  the  other  one  to  the  site,  which  were  designed  to  approximately  equal¬ 
ize  all  peak  amplitudes.  This  method  is  in  agreement  with  the  multiplicative  structure  of  the 
spectra,  and  the  exact  values  of  the  equalizing  factors  do  not  matter.  The  equalizing  of  traces 
conflicts  some  with  the  optimization  of  the  S/N  ratios  in  the  output,  but  as  mentioned  above 
this  is  not  our  main  goal.  After  processing,  the  initial  normalization  factors  can  be  divided  out 
to  estimate  the  proper  relative  amplitudes  of  the  deconvolved  traces. 

The  iteration  procedure  is  started  with  the  initialization  of  the  time  domain  representations 
of  the  source  and  site  factors  to  delta  functions.  The  source  factors  are  estimated  first  by  sum¬ 
ming  the  deconvolved  results  over  the  sites.  This  way  the  first  estimates  of  source  functions 
are  similar  to  those  obtained  by  Marshall  (1972),  but  during  the  subsequent  iterations  both  the 
source  and  site  factors  have  been  further  refined  by  factoring. 

In  order  to  ensure  the  desired  deconvolution  effect  we  removed  the  division  by  the  noise 
spectra,  which  would  have  resulted,  especially  for  NTS  events  which  have  the  maxima  in  S/N 
ratio  below  1  Hz,  in  no  visible  deconvolution  at  all.  Instead,  we  have  assigned  a  unit  weight 
to  all  frequencies  which  had  S/N  ratios  above  1  and  excluded  the  frequency  ranges  with  poor 
S/N  ratios  from  the  iterations.  Again,  this  degrades  the  S/N  in  the  output,  but  we  accept  this 
as  the  price  of  increased  resolution. 
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Table  1 


Source  Parameters  for  NTS  Events  Recorded  at  NORSAR 

Date 

Name 

Location 

mb  (ISC) 

14  May  1975 

Tybo 

Pahute  Mesa 

5.9 

03  Jun  1975 

Stilton 

Pahute  Mesa 

5.8 

19  Jun  1975 

Mast 

Pahute  Mesa 

5.9 

28  Oct  1975 

Kasseri 

Pahute  Mesa 

6.2 

20  Nov  1975 

Inlet 

Pahute  Mesa 

5.9 

Table  3 


Source  Parameters  for  Eastern  Kazakh  Events 
Recorded  at  EKA 

Date 

Location 

mb  (ISC) 

15  Jan  1965 

Eastern  KTS 

6.3 

27  Apr  1975 

Eastern  KTS 

5.6 

4  Jul  1976 

Eastern  KTS 

5.8 

7  Dec  1976 

Eastern  KTS 

5.9 

11  Jun  1978 

Eastern  KTS 

5.9 

KTS=Eastem  Kazakh  Test  Site 

Eastern  KTS  corresponds  to  a  testing  sites  within  the  KTS  as  identified  by  Dahlman  and 
Israelson  (1977,  pp  182-183). 


Table  4 


Source  Parameters  for  Astrakhan  Events 
Recorded  at  the  RSTN 


Date 

Origin  Time 

"H,  (ISC) 

6  Oct  1982 

05:59:57.4 

5.2 

6  Oct  1982 

06:04:57.4 

5.2 

6  Oct  1982 

06:09:57.4 

5.2 

6  Oct  1982 

06:14:57.4 

5.4 

Parameters  used  in  Multichannel  Deconvolutions 


Source 

Receiver 

Number  of 

Number  of 

Frequency 

Sample 

Region 

Array 

Sources 

Receivers 

Bandwidth 

t* 

Rate 

(Hz) 

(sec) 

(sec'1) 

NTS 

NORSAR 

5 

5 

0.4-3.5 

0.45 

20 

KTS 

NORSAR 

8 

7 

0.4-5.0 

0.15 

20 

KTS 

EKA 

5 

12 

0.4-5.0 

0.15 

20 

Astrakhan 

RSTN 

4 

4 

0. 6-5.0 

0.15 

40 
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Figures  11  to  14  show  the  deconvolved  source-time  functions  for  the  four  data  sets. 
Most  of  the  buried  explosions  show  phases  which  could  readily  be  interpreted  as  pP,  though  all 
the  waveforms  show  a  good  deal  more  complexity  than  can  be  described  simply  by  P  +  pP.  A 
good  deal  of  energy  arrives  for  a  couple  of  seconds  after  the  "P  +  pP"  in  nearly  all  the  events. 
Average  estimates  of  pP  delay  times  are  0.7  sec  for  NTS  events  recorded  at  NORSAR  (Figure 
11),  0.55  sec  for  KTS  events  recorded  at  NORSAR  (Figure  12),  0.4  sec  for  KTS  events 
recorded  at  EKA  (Figure  13),  and  1.0  sec  for  Astrakhan  events  recorded  at  RSTN  stations 
(Figure  14).  In  general,  the  KTS  events  (Figures  12  and  13)  have  simpler  source  time  func¬ 
tions  than  the  NTS  events  (Figure  11).  From  figures  12  and  13,  there  is  no  obvious  difference 
in  the  pP  delay  times  between  shots  at  the  eastern  and  central  test  sites  of  the  KTS  (see  pages 
182  and  183  of  Dahlman  and  Israelson  [1977]  for  a  description  of  the  divison  of  the  KTS  into 
eastern,  central,  and  western  testing  sites).  The  presumed  cratering  event  of  15  January  1965 
(bottom  of  Figure  13)  is  a  notable  exception  to  the  other  events.  The  source  waveform  is  visi¬ 
bly  more  complex  than  those  associated  with  buried  explosions  at  the  eastern  KTS  and  the 
definition  of  pP  is  much  more  ambiguous.  The  Astrakhan  events  (Figure  14)  also  have  fairly 
complex  waveforms  (perhaps  due  to  the  sharply  contrasting  velocity  structure  in  a  salt  dome 
environment)  and  long  inferred  pP  times,  consistent  with  suggestions  of  deep  burial. 

There  is  remarkable  agreement  between  the  deconvolved  source  terms  from  these  decon¬ 
volutions  and  from  the  stacked  deconvolutions  presented  in  the  previous  section.  In  particular, 
compare  Figure  7  with  Figure  11  and  Figure  8  with  Figure  13. 

The  deconvolved  source  terms  for  the  KTS  explosions  at  EKA  (Figure  13)  were  used  to 
estimate  the  magnitude  bias  between  the  four  buried  explosions  and  the  cratering  shot.  In  an 
effort  to  compensate  for  the  differences  in  yields  between  the  five  events,  the  deconvolved 
source  terms  were  normalized  so  that  the  first  zero  to  peak  swing  of  the  P  wave  has  an  ampli- 
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Figure  11.  Deconvolved  source  terms  for  NTS  events  recorded  at  NORSAR.  Each  trace  is 
labeled  with  the  event  name  and  the  locations  of  the  P  and  pP  (if  observed)  phases.  Event 
parameters  for  these  events  are  given  in  Table  1. 
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Figure  12.  Deconvolved  source  terms  for  Eastern  Kazakh  Test  Site  (KTS)  events  recorded  at 
NORSAR.  Each  trace  is  labeled  with  the  event  name  and  the  locations  of  the  P  and  pP  (if 
observed)  phases.  Event  parameters  for  these  events  are  given  in  Table  2. 
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Figure  13.  Deconvolved  source  terms  for  KTS  events  recorded  at  EKA.  Each  trace  is  labeled 
with  the  event  name  and  the  locations  of  the  P  and  pP  (if  observed)  phases.  Event  parameters 
for  these  events  are  given  in  Table  3.  The  bottom  trace  is  a  presumed  cratering  event 
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Figure  14.  Deconvolved  source  terms  for  Astrakhan  events  recorded  at  the  RSTN.  Each  trace 
is  labeled  with  the  event  name  and  the  locations  of  the  P  and  pP  (if  observed)  phases.  Event 
parameters  for  these  events  are  given  in  Table  4. 
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tude  of  1.  Then  each  trace  was  convolved  with  a  wavelet  derived  from  a  von  Seggem  and 
Blandford  source  function,  t*  =  0.15  sec,  and  an  EKA  instrument  response.  An  example  of 
such  filtered  source  terms  is  shown  in  Figure  15  for  a  yield  of  125  kt.  Thus,  the  magnitude 
bias  is  estimated  from  the  differences  in  waveforms  due  to  the  differences  in  the  deconvolved 
source  time  functions  between  the  cratering  and  non-cratering  events.  The  magnitude  bias  was 
calculated  for  three  yields,  60,  125,  and  300  kt,  in  two  different  ways: 


Amb  =  mb(n— cr)  -  mb(cr)  =  log 


*c,cr 


(6) 


*  « 
Ac 

•  * 
Ac 

Amb  =  log 

aT 

n-cr 

log 

a7 

where  Aa  and  Ac  are  the  "a"  phase  and  "c”  phase  amplitudes,  respectively,  and  cr  and  n-cr 
refer  to  cratering  and  non-cratering  events,  respectively.  The  results  are  shown  in  Figure  16, 
and  Amb(n-cr  -  cr)  -  0.09  to  0.22  magnitude  units.  It  is  interesting  that  the  results  from  using 
Equation  (7)  appear  more  robust  than  those  from  using  Equation  (6).  This  may  be  because  the 
yield  differences  for  the  events,  which  were  only  partly  corrected  for  by  normalizing  each 
trace,  are  in  a  sense  "factored  out"  when  the  ratios  of  the  "a"  and  "c"  phase  amplitudes  are 
used. 


The  deconvolved  site  terms  are  a  by  product  of  the  deconvolution  process  in  the  sense 
that  we  are  mostly  interested  in  the  source  time  functions  in  this  study.  The  site  terms  are 
shown  in  Figures  B1  to  B4  of  Appendix  B;  they  are  "wrapped  around"  because  the  deconvolu¬ 
tion  uses  finite  time  series.  There  are  substantial  variations  in  the  site  terms,  even  among  sta¬ 
tions  near  to  each  other  (such  as  those  in  the  NORSAR  and  EKA  arrays),  and  all  of  the  site 
terms  appear  to  be  the  source  of  some  of  the  P  wave  coda.  To  show  the  effect  of  the  site 
terms  on  seismic  waves  in  more  detail,  Figures  17  to  19  show  site  terms  that  have  been 
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Figure  15.  Deconvolved  source  terms  from  the  set  of  eastern  KTS  data  recorded  at  EKA 
which  have  been  filtered  with  a  wavelet  consisting  of  a  von  Seggem  and  Blandford  pulse  of 
yield  125  kt,  a  t*  =  0.15  sec,  and  an  EKA  instrument  response.  The  bottom  trace  is  from  the 
eastern  KTS  cratering  shot. 
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Figure  16.  Magnitude  difference  estimates  between  the  presumed  cratering  event  at  eastern 
KTS  and  four  buried  eastern  KTS  explosions  for  yields  of  60,  125,  and  300  kL  The  magni¬ 
tude  differences  were  calculated  from  equations  (6)  and  (7)  in  the  text.  Aa  and  Ac  are  the  "a" 
phase  and  "c"  phase  amplitudes,  respectively,  and  cr  and  n-cr  refer  to  cratering  and  non¬ 
cratering  events,  respectively. 


Figure  17.  Examples  of  filtered  site  terms  for  the  set  of  NTS  data  recorded  at  NORSAR.  In 
each  pair,  the  top  trace  is  the  deconvolved  site  term,  which  has  been  filtered  in  the  bottom 
trace  with  a  von  Seggern  and  Blandford  wavelet  assuming  a  yield  of  500  kt,  t*  =  0.45  sec,  and 
a  NORSAR  instrument  response. 
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Figure  18.  Examples  of  filtered  site  terms  for  the  set  of  KTS  data  recorded  at  EKA.  In  each 
pair,  the  top  trace  is  the  deconvolved  site  term,  which  has  been  filtered  in  the  bottom  trace  with 
a  von  Seggern  and  Blandford  wavelet  assuming  a  yield  of  130  kt,  t*  =  0.15  sec,  and  an  EKA 
instrument  response. 
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convolved  with  von  Seggem  and  Blandford  wavelets  including  t*  and  an  appropriate  instru¬ 
ment  response.  In  Figure  19,  the  predicted  arrival  time  of  PcP  from  Astrakhan  has  been 
marked.  The  influence  of  the  site  function  on  the  coda  is  clearly  seen,  as  is  the  difference 
between  sites  that  are  closely  located  geographically.  Also,  three  NORSAR  sites  were  used  in 
studies  of  data  from  events  at  two  different  azimuths  (NTS  and  KTS).  The  site  terms  for  these 
receivers  from  the  two  azimuths  are  shown  in  Figure  20.  The  initial  response  of  each  site 
(unction  is  similar  for  the  events  from  different  azimuths,  while  there  is  more  variation  in  the 
coda,  especially  in  NB5  (bottom  of  Figure  20). 

The  resolution  kernels  for  the  deconvolutions  are  shown  in  Figures  B5  to  B8  of  Appen¬ 
dix  B.  The  width  of  the  resolution  kernels  depends  on  the  frequency  bandwidth  of  the  decon¬ 
volutions;  for  the  source-receiver  pairs  studied  in  this  report,  the  width  is  about  0.15  sec  at  half 
maximum. 

As  one  check  on  the  results,  reconstructions  were  calculated  by  convolving  the  appropri¬ 
ate  source  wavelet,  deconvolved  source  time  function,  and  deconvolved  site  time  function.  All 
of  the  original  traces  are  in  Appendix  C,  shown  for  comparison  with  the  reconstructions  in 
Appendix  D.  Figures  21  and  22  show  examples  of  traces  from  an  NTS  event  and  a  KTS  event 
along  with  their  reconstructions.  The  similarity  of  the  observed  and  reconstructed  traces  is 
very  good,  better  in  fact  than  some  reported  reconstructions  made  with  "relative  receiver  func¬ 
tions",  probably  because  we  use  an  array  of  small  areal  extent  rather  than  widely  spaced  sta¬ 
tions.  Moreover,  the  correlation  coefficients  between  the  original  and  reconstructed  traces  are 
commonly  above  0.9  for  the  whole  25.6  second  length  of  the  samples  chosen  (12.8  sec  for  the 
Astrakhan  data).  Also  note  how  well  the  details  of  the  coda  have  been  preserved  in  the  recon¬ 
structions.  The  reconstructions  prove,  at  least,  that  the  observed  waveforms  are  not  incon¬ 
sistent  with  some  plausible  pP  parameters  coupled  with  low  t*  that,  unlike  values  used  in  pre- 
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Figure  21.  Five  traces  from  the  NTS  explosion  Mast  as  recorded  at  NORSAR  (top  of  each 
pair  of  traces),  along  with  their  reconstructions  (bottom  of  each  pair  of  traces).  The  receiver 
names  are  given  to  the  left  of  each  pair  of  traces,  as  is  the  correlation  coefficient  between  each 
original  trace  and  its  reconstruction. 
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Figure  22.  Five  traces  from  the  eastern  KTS  cratering  shot  of  15  January  1965  as  recorded  at 
EKA  (top  of  each  pair  of  traces),  along  with  their  reconstructions  (bottom  of  each  pair  of 
traces).  The  receiver  names  are  given  to  the  left  of  each  pair  of  traces,  as  is  the  correlation 
coefficient  between  each  original  trace  and  its  reconstruction. 
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vious  synthetic  simulations,  also  agree  with  spectral  data. 


Conclusions 


Most  of  the  ambiguities  associated  with  the  pP  estimation  methods  proposed  in  the  past 
are  associated  with  the  fact  that,  for  most  explosions,  the  simple  P+pP  model  is  not  valid. 
This  can  be  verified  by  analyses  of  interevent  spectral  ratios  and  the  results  of  our  multichannel 
deconvolutions.  In  addition,  for  single  channels  the  site  effects  will  also  disguise  any  spectral 
modulation  due  to  P-pP  interference  even  if  the  P+pP  model  is  valid.  Finite  difference  simula¬ 
tions  of  P  waves  from  explosion  sources  emplaced  in  laterally  heterogeneous  structures  also 
support  the  idea  that  complex  arrival  patterns,  rather  than  simple  P+pP  sequences,  should  be 
common. 

Multichannel  deconvolutions  of  four  data  sets  consisting  of  array  recordings  of  suites  of 
events  at  several  test  sites  gave  the  following  results: 

*  NTS  events  deconvolved  at  NORSAR  show,  in  general  complex  source  functions 
not  easily  interpretable  in  terms  of  the  P+pP  model. 

*  Eastern  Kazakh  explosions  deconvolved  at  EKA  and  NORSAR  appear  to  have 
simpler  source  functions  with  an  identifiable  pP  arrival  with  a  time  delay  in  the  0.3- 
0.6  sec  range. 

*  The  deconvolved  source  function  of  a  presumed  cratering  explosion  at  Eastern 
Kazakh  is  visibly  different  from  the  deconvolved  source  functions  for  buried  explo¬ 
sions  at  the  same  test  site.  Assuming  the  first  positive  spike  to  be  direct  P,  the 
effect  of  the  cratering  on  mb  was  estimated  as  a  reduction  of  0.09  to  0.22  magnitude 
units. 

*  For  a  set  of  Astrakhan  events  deconvolved  using  the  RSTN  stations  we  see  a  late 
arriving  pP  1  sec  after  the  direct  P,  as  well  as  other  late  arrivals. 

*  Reconstructions  from  the  site  and  source  time  series  and  the  source  wavelets  show 
excellent  agreement  with  the  original  traces  including  fine  details  in  the  coda.  The 
coda  is  partially  due  to  both  near  site  and  source  scattering. 
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These  results  have  profound  implications  to  the  study  of  scattering  of  seismic  waves  near 
the  sources  and  receivers.  Since  most,  but  not  all,  of  the  source  factors  are  associated  with 
relatively  simple  time  series  compared  to  the  site  factors,  the  associated  events  must  have  gen¬ 
erated  little  scattered  energy  near  the  source.  If  this  were  not  true  we  would  not  have  been 
able  to  reconstruct  the  associated  time  series.  We  cannot  think  of  any  conceivable  physical 
mechanism  that  would  generate  complex  and  similar  source  wavetrains  for  explosions  at  vari¬ 
ous  locations  within  the  same  test  site. 

Suggested  Future  Work 

Although  in  the  work  presented  we  attempted  to  ensure  the  uniqueness  of  our  results  by 
using  reasonable  starting  assumptions,  much  more  can  be  done  to  reduce  the  uncertainties  in 
the  final  results.  There  are  several  possible  constraints,  all  based  on  common  sense  arguments 
that  have  not  been  as  yet  used  in  our  processing  of  the  data.  For  instance,  one  might  use  the 
expectation  that  the  time  domain  representations  of  site  factors  be  uncorrelated  beyond  the  ini¬ 
tial  few  seconds.  Similarly  one  might  constrain  the  later  parts  of  the  source  time  functions  to 
be  uncorrelated.  One  might  also  exploit  the  fact  that  some  sites  show  waveforms  that  are  more 
coherent  with  each  other,  shorter  and  simpler  than  at  some  "anomalous"  sites.  This  does  not 
mean  however  that,  similar  to  the  "transparent"  station  concept,  there  is  no  site  distortion  at 
these  sites.  Rather  it  may  simply  mean  that  the  effects  of  the  scattering  operator  are  similar  in 
the  first  second  and  less  severe  at  later  times.  We  must  also  point  out  that  the  iterative,  fre¬ 
quency  domain  algorithm  for  this  process  is  by  no  means  essential  and  other  algorithms  may 
be  just  as  viable.  Successive  time  domain  deconvolutions  with  one-sided  causal  filters  may  also 
be  performed  iteratively,  although  such  filters  may  be  less  effective  in  removing  complex  site 
effects.  On  the  other  hand,  causal  outputs  are  more  appealing  and  appear  more  realistic. 
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Despite  the  fact  that  the  process  we  have  applied  is  not  constrained  to  be  causal,  the  source 
and  site  time  series  as  well  as  the  reconstructions  do  not  show  appreciable  precursors.  A  better 
understanding  of  the  convergence  properties  of  the  iterative  method  would  also  be  desirable. 
We  prefer  the  iterative  approach  because  it  allows  us  to  change  constraints  as  we  go  along,  but 
the  idea  of  factoring  can  be  exploited  without  any  iteration.  Most  of  the  effort  in  this  study 
was  expanded  in  developing,  implementing  and  testing  the  algorithm  and  little  time  was  left  for 
"fine  tuning"  and  optimization.  The  results  shown  here  are  probably  inferior  to  the  best  that 
can  be  produced  by  the  method  after  optimizing  the  weighting  of  frequencies  and  the  initial 
conditions. 

After  testing  out  and  optimizing  the  method,  the  implementation  of  an  intelligent,  adap¬ 
tive  scheme  for  routine  determination  of  source  time  functions  at  seismic  arrays  becomes  a 
practical  possibility.  The  deconvolution  process  also  easily  lends  itself  for  incorporation  into 
"intelligent"  learning  systems  which  adapt  themselves  as  new  data  are  added.  For  instance,  if 
we  are  satisfied  with  the  previous  estimates  of  the  site  factors  we  can  "freeze"  these  and  the 
process  can  be  used  for  routine  estimation  of  the  source  time  functions  from  the  same  test  site. 
If  we  have,  from  other  sources  or  methods,  "true”  values  of  pP  parameters  then  we  can  keep 
the  associated  source  function  constant  and  rerun  the  process  to  obtain  estimates  of  site  and 
source  factors  conditional  to  such  constraints.  Improvements  in  t*,  explosion  source  models 
etc.  can  be  naturally  incorporated  in  the  scheme,  and  the  results  produced  in  this  report  must 
thus  not  be  considered  as  to  be  the  best  that  can  be  obtained.  The  algorithm  could  be  installed 
permanently  in  a  computer  and  all  data  originating  from  selected  test  sites  may  be  processed. 
Instead  of  retaining  all  the  past  data  only  the  updated  site  factors  would  be  kept  for  the  pro¬ 
cessing  of  new  signals  along  various  azimuths.  These  may  take  the  place  of  stored  beam 
parameters.  It  would  be  possible,  using  multiple  arrays,  to  look  for  effects  of  strain  release 
along  various  azimuths  and  similarly  to  the  comparison  of  cratering  and  noncratering  events  in 
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Appendix  A 

Spectral  ratios  of  pairs  of  NTS  events  at  NORSAR  stations. 
Spectral  ratios  of  NTS  events  at  pairs  of  NORSAR  stations. 


Table  A1 


Parameters  used  in  Forming 

Synthetic  Spectral  Ratios 

Event 

pP-P  (sec) 

pP/P  Yield  (kt) 

Camembert  1.15 

0.9 

1003 

Inlet 

0.95 

1.0 

315 

Kasseri 

1.05 

0.9 

1272 

Mast 

0.9 

0.9 

579 

Muenster 

1.25 

0.8 

1133 

Stilton 

0.85 

0.9 

222 

Tybo 

0.9 

1.3 

386 

From  Lay  et  al  (1984a).  pP-P  and  pP/P  are  from  Mast  as  the  master  event;  the  yields  are 
the  average  of  the  yields  from  the  four  different  master  events. 
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Figure  Al.  Spectral  ratios  between  NTS  events  Inlet  and  Kasseri  at  NORSAR  stations. 


53 


FREQUENCY  (Hz) 

Figure  A2.  Spectral  ratios  between  NTS  events  Inlet  and  Mast  at  NORSAR  stations. 
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Appendix  B 


Deconvolved  site  time  functions. 
Resolution  kernels  from  deconvolutions. 
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Figure  Bl.  Deconvolved  site  terms  for  NTS  events  recorded  at  NORSAR. 
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Figure  B2.  Deconvolved  site  terms  for  KTS  events  recorded  at  NORSAR. 
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Figure  B2  (continued). 
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Figure  B3.  Deconvolved  site  terms  for  KTS  events  recorded  at  EKA. 
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Figure  B3  (continued). 


E-AD 


E-AD 


A  o  0  S  E  L 


SITE  FUNCTION* 


Figure  B3  (continued). 
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Figure  B5.  Resolution  kernels  for  NTS  events  recorded  at  NORSAR. 
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Figure  B6  (continued). 
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Figure  B7.  Resolution  kernels  for  KTS  events  recorded  at  EKA. 


(THIS  PAGE  INTENTIONALLY  LEFT  BLANK) 


>WH 


m 


»  .  fo  lf.  <  ■ 


75 


'  VV  V./ ./> 

vl’.r'vi  VN.“^ 


y.v. 


,*  V  „•  '.“  V  V 


v:a'-.  •■ 

V  -j 


a 


'  V 


KASSERI 


81 


I 


«! 


*nO  SEC 


22MAR71 


84 


*  '  , 


■f  r  «- 


I02C 


J07B 


SEC  30JUN71 


86 


\  /,v. 
»  — 


t 


►  /, 


XX 


16FEB73 


Hip  VyijWVy  V Yh^V^^Atv^TV ^-^AjxNvv'' vwv^V^aJ 


■,/A^/y\/'yWA'VKrT\/V^^yyWvAA/V\yvA/V\V"^v^/^^'^ 


HoO  SEC 


7 AUG 7 5 


97 


N02C 


N07B 


J  VI 


H.O  SEC 


7AUG75 


98 


V^/WvV 


WV^\/VV^v -^/^V ' 


l/VV 


'A^V'v/s/  V~ 


Ho  0  SEC 


7 DEC 76 


Ho  0  SEC 


7DEC76 


f|/\/\/y\A/v^\^^/VA- 


1  1  JUN78 


'■'v jJ 


wtnAhfJi 


Appendix  D 

Reconstructed  traces  of  explosion  data. 
NTS  data  recorded  at  NORSAR. 
Kazakh  data  recorded  at  NORSAR. 
Kazakh  data  recorded  at  EKA. 
Azgir  data  recorded  at  the  RSTN. 
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